/* Computer Graphics, Assignment, Volume rendering with cubes/points/isosurface
 *
 * Student name .... Dexter Drupsteen
 * Student email ... dexter.drupsteen@gmail.com
 * Collegekaart .... 5887224
 * Student name .... Jari Ferguson
 * Student email ... jari@jsferguson.nl
 * Collegekaart .... (UvAID)6375847
 * Date ............ 23/11/2011
 *
 * (always fill in these fields before submitting!!)
 */

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include "marching_tetrahedra.h"

/* Compute a linearly interpolated position where an isosurface cuts
   an edge between two vertices (p1 and p2), each with their own
   scalar value (v1 and v2) */

static vec3
interpolate_points(unsigned char isovalue, vec3 p1, vec3 p2, unsigned char v1, unsigned char v2)
{
    /* Initially, simply return the midpoint between p1 and p2.
       So no real interpolation is done yet */

	vec3 vec_tmp = v3_add(v3_multiply(p1,(abs(isovalue - v1) / abs(v2 - v1))), v3_multiply(p2, 1 - (abs(isovalue - v1) / abs(v2 - v1))));
    return v3_create(vec_tmp.x * sizex, vec_tmp.y * sizey, vec_tmp.z * sizez);
}

/* Using the given iso-value generate triangles for the tetrahedron
   defined by corner vertices v0, v1, v2, v3 of cell c.

   Store the resulting triangles in the "triangles" array.

   Return the number of triangles created (either 0, 1, or 2).

   Note: the output array "triangles" should have space for at least
         2 triangles.
*/

static int
generate_tetrahedron_triangles(triangle *triangles, unsigned char isovalue, cell c, int v0, int v1, int v2, int v3)
{
	//associate with vertex
	unsigned char vert_bit = 0;
	int nr_triangles = 0;
	if (c.value[v0] > isovalue )
	{
		vert_bit = 1;
	}
	else vert_bit = 0;
	
	if (c.value[v1] > isovalue )
	{
		vert_bit = vert_bit + 2;
	}
	
	if (c.value[v2] > isovalue )
	{
		vert_bit = vert_bit + 4;
	}
	
	if (c.value[v3] > isovalue )
	{
		vert_bit = vert_bit + 8;
	}
	
	// start generating the triangles
	if (vert_bit == 1 || vert_bit == 14)
	{
		triangles[0].p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], c.value[v0], c.value[v1]);
		triangles[0].p[1] = interpolate_points(isovalue, c.p[v0], c.p[v3], c.value[v0], c.value[v3]);
		triangles[0].p[2] = interpolate_points(isovalue, c.p[v0], c.p[v2], c.value[v0], c.value[v2]);
		nr_triangles = 1;
	}
	
	if (vert_bit == 2 || vert_bit == 13)
	{
		triangles[0].p[0] = interpolate_points(isovalue, c.p[v1], c.p[v0], c.value[v1], c.value[v0]);
		triangles[0].p[1] = interpolate_points(isovalue, c.p[v1], c.p[v3], c.value[v1], c.value[v3]);
		triangles[0].p[2] = interpolate_points(isovalue, c.p[v1], c.p[v2], c.value[v1], c.value[v2]);
		nr_triangles = 1;
	}
	
	if (vert_bit == 4 || vert_bit == 11)
	{
		triangles[0].p[0] = interpolate_points(isovalue, c.p[v2], c.p[v0], c.value[v2], c.value[v0]);
		triangles[0].p[1] = interpolate_points(isovalue, c.p[v2], c.p[v3], c.value[v2], c.value[v3]);
		triangles[0].p[2] = interpolate_points(isovalue, c.p[v2], c.p[v1], c.value[v2], c.value[v1]);
		nr_triangles = 1;
	}
	
	if (vert_bit == 8 || vert_bit == 7)
	{
		triangles[0].p[0] = interpolate_points(isovalue, c.p[v3], c.p[v0], c.value[v3], c.value[v0]);
		triangles[0].p[1] = interpolate_points(isovalue, c.p[v3], c.p[v2], c.value[v3], c.value[v2]);
		triangles[0].p[2] = interpolate_points(isovalue, c.p[v3], c.p[v1], c.value[v3], c.value[v1]);
		nr_triangles = 1;
	}
	
	if (vert_bit == 3 || vert_bit == 12)
	{
		triangles[0].p[0] = interpolate_points(isovalue, c.p[v0], c.p[v3], c.value[v0], c.value[v3]);
		triangles[0].p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], c.value[v0], c.value[v2]);
		triangles[0].p[2] = interpolate_points(isovalue, c.p[v1], c.p[v2], c.value[v1], c.value[v3]);
		triangles[1].p[0] = interpolate_points(isovalue, c.p[v1], c.p[v3], c.value[v1], c.value[v3]);
		triangles[1].p[1] = interpolate_points(isovalue, c.p[v1], c.p[v3], c.value[v1], c.value[v2]);
		triangles[1].p[2] = interpolate_points(isovalue, c.p[v0], c.p[v2], c.value[v0], c.value[v2]);
		nr_triangles = 2;
	}
	
	if (vert_bit == 5 || vert_bit == 10)
	{
		triangles[0].p[0] = interpolate_points(isovalue, c.p[v3], c.p[v0], c.value[v3], c.value[v0]);
		triangles[0].p[1] = interpolate_points(isovalue, c.p[v3], c.p[v2], c.value[v3], c.value[v2]);
		triangles[0].p[2] = interpolate_points(isovalue, c.p[v1], c.p[v2], c.value[v1], c.value[v2]);
		triangles[1].p[0] = interpolate_points(isovalue, c.p[v1], c.p[v0], c.value[v1], c.value[v0]);
		triangles[1].p[1] = interpolate_points(isovalue, c.p[v1], c.p[v2], c.value[v1], c.value[v2]);
		triangles[1].p[2] = interpolate_points(isovalue, c.p[v0], c.p[v3], c.value[v0], c.value[v3]);
		nr_triangles = 2;
	}
	
	if (vert_bit == 6 || vert_bit == 9)
	{
		triangles[0].p[0] = interpolate_points(isovalue, c.p[v3], c.p[v2], c.value[v3], c.value[v2]);
		triangles[0].p[1] = interpolate_points(isovalue, c.p[v3], c.p[v1], c.value[v3], c.value[v1]);
		triangles[0].p[2] = interpolate_points(isovalue, c.p[v0], c.p[v2], c.value[v0], c.value[v2]);
		triangles[1].p[0] = interpolate_points(isovalue, c.p[v1], c.p[v3], c.value[v1], c.value[v3]);
		triangles[1].p[1] = interpolate_points(isovalue, c.p[v1], c.p[v0], c.value[v1], c.value[v0]);
		triangles[1].p[2] = interpolate_points(isovalue, c.p[v0], c.p[v2], c.value[v0], c.value[v2]);
		nr_triangles = 2;
	}
	
	if (nr_triangles == 1)
	{
		vec3 vector1, vector2;
		
		for (int i = 0; i < 3; i++)
		{
			vector1 = v3_subtract(triangles[0].p[(i + 1) % 3], triangles[0].p[i]);
			vector2 = v3_subtract(triangles[0].p[(i + 2) % 3], triangles[0].p[i]);
			triangles[0].n[i] = v3_normalize(v3_crossprod(vector1, vector2));	
		}
		
	}
	if (nr_triangles == 2)
	{
		vec3 vector1, vector2, vector3, vector4;
		
		for (int i = 0; i < 3; i++)
		{
			vector1 = v3_subtract(triangles[0].p[(i + 1) % 3], triangles[0].p[i]);
			vector2 = v3_subtract(triangles[0].p[(i + 2) % 3], triangles[0].p[i]);
			triangles[0].n[i] = v3_normalize(v3_crossprod(vector1, vector2));
			
			vector3 = v3_subtract(triangles[1].p[(i + 1) % 3], triangles[1].p[i]);
			vector4 = v3_subtract(triangles[1].p[(i + 2) % 3], triangles[1].p[i]);
			triangles[1].n[i] = v3_normalize(v3_crossprod(vector3, vector4));	
		}
	}
	
	return nr_triangles;
    
}

/* Generate triangles for a single cell for the given iso-value. This function
   should produce at most 6 * 2 triangles (for which the "triangles" array should
   have enough space).

   Use calls to generate_tetrahedron_triangles().

   Return the number of triangles produced
*/

int
generate_cell_triangles(triangle *triangles, cell c, unsigned char isovalue)
{
    int nr_produced_triangles = 0;
    int length = 0;
    
    triangle new_triangles[2];
    
    length = generate_tetrahedron_triangles(new_triangles, isovalue, c, 0,1,3,7);
    for (int i = 0; i < length; i++)
    {
    	triangles[nr_produced_triangles++] = new_triangles[i];
    }
    
    length = generate_tetrahedron_triangles(new_triangles, isovalue, c, 0,2,6,7);
    for (int i = 0; i < length; i++)
    {
    	triangles[nr_produced_triangles++] = new_triangles[i];
    }
    
    length = generate_tetrahedron_triangles(new_triangles, isovalue, c, 0,1,5,7);
    for (int i = 0; i < length; i++)
    {
    	triangles[nr_produced_triangles++] = new_triangles[i];
    }
    
    length = generate_tetrahedron_triangles(new_triangles, isovalue, c, 0,2,3,7);
    for (int i = 0; i < length; i++)
    {
    	triangles[nr_produced_triangles++] = new_triangles[i];
    }
    
    length = generate_tetrahedron_triangles(new_triangles, isovalue, c, 0,4,5,7);
    for (int i = 0; i < length; i++)
    {
    	triangles[nr_produced_triangles++] = new_triangles[i];
    }
    
    length = generate_tetrahedron_triangles(new_triangles, isovalue, c, 0,4,6,7);
    for (int i = 0; i < length; i++)
    {
    	triangles[nr_produced_triangles++] = new_triangles[i];
    }
    return nr_produced_triangles;
}
